Airy-like pulses in models of large molecular chains, and 
conservative numerical methods for quasi-linear Hamiltonian 

systems 

Brenton LeMesurier 
Department of Mathematics, College of Charleston, South Caroling 

(Dated: January 15, 2013) 

Abstract 

The phenomenon of coherent energetic pulse propagation in macromolecular chains such as a- 
helix protein is studied using the Davydov-Scott model, with both numerical studies using a new 
unconditionally stable fourth order accurate energy-momentum conserving time discretization, and 
with analysis based on ideas of center manifold theory. 

It is shown that for physically natural impulsive initial data, the coherent traveling pulses seen 
have a form related to the Airy function, but with rapid variation of phase along the chain. This 
can be explained in terms of a new continuum limit approximation by the third derivative nonlinear 
Schrodinger equation, which differs from the previous continuum limit approximations related to 
the standard NLS equation. 

A theorem is given describing the construction of such conservative time discretizations for a 
large class of Hamiltonian systems. 
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I. INTRODUCTION 



A. Exciton-Oscillator Models of Macromolecular Chains 

Various phenomena in macromolecular chains, such as proteins and DNA, are modeled 
by large systems of ODE's in Hamiltonian form, with conservation of energy and possibly 
some quadratic invariants, here called "momenta" . In particular, many such models, along 
with spatial discretizations of various nonlinear wave equations, involve exciton- oscillator 
systems 



S V 



(1) 



V 



with symmetry conditions L_s = Lg, Vl^s = ^s- These have Hamiltonian form 

i 
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dq^ dH dpy dH 



(4) 



dt dpy ' dt dq^ 
Note that ipn and ip^ are here treated formally as independent variables, but in practice 

they are complex conjugates, the Hamiltonian is real, and so only the first equation in ([3]) 

is needed. 

There is also a related class of Lattice Nonlinear Schrodinger [LNLS] equation models 

+ ^^'^^ri+s = 2^Xs\^n+s\'^^n (5) 
s s 

with Hamiltonian 

'H = ^Lsij^llJn+s + ^Xs\^n\'^\^n+s\'^- (6) 

n,s n.s 

Section |IT1 introduces the main example considered here, one developed by Alexandr 
Davydov, Alwyn Scott et al [lllsl for energy propagation in a-helix protein, which is the 
dominant helical form occurring in parts of protein molecules 



dt 

dt^ 



Jqn 

fn-j-r- - {q„-3 - 2g„ + qn+3) = l^nP - IV'n-sT) 



along with a lattice NLS approximation thereof. 
.dip. 



dt 



+ J(^/'„+3 + ^n-s) - L{^n+1 + ^n-l) + = 0. 



(8) 



The physical meaning and parameter values are discussed further in the next section. 

Many of the methods and results here also apply to related models as described in surveys 
by Davydov {g] and Scott j?], and considered for example in the more recent work of Ivic et 



al m 



10|. 



B. Numerical Observations and Analysis 

Section IIIII introduces a new fourth order accurate time-discretization method that ex- 
actly conserves the Hamiltonian ("energy") and all quadratic invariants ("momenta"), with 
unconditional stability to handle the stiffness that can arise. A difference calculus proce- 
dure is described for constructing such time discretizations with conservation of all quadratic 
invariants as well as the energy (Hamiltonian). 

This discrete gradient method adapts easily to model refinements such as stochastic 
modeling of interactions with the aqueous environment and more accurate modeling of intra- 
molecular forces, and to related questions such as signal propagation in DNA. The order of 
accuracy can easily be increased when this is desirable. 

Section IIVI then applies this method in numerical studies that partially confirm previous 
observations and conjectures that sustained pulses can carry energy further and with less 
dispersion than one might expect from linearized approximations. However, a previous 
continuum limit approximation in terms of the nonlinear Schrodinger [NLS] equation and 
its sech soliton traveling wave solutions is shown to be invalid for typical initial data, due 
to a fast variation of the phase along the chain. 

Instead, it is shown in Section |V] that the coherent pulse phenomenon seen is in part ap- 
proximated by a new continuum limit, in terms of the third derivative nonlinear Schrodinger 
equation, having the Airy PDE as its linear part. The pulse structure is related to solutions 
of it'd linearization, the Airy PDE, which are in terms of the Airy function. With larger 
initial data, and thus stronger nonlinear effects, the pulse becomes more unimodal, with 
great concentration of signal strength in leading pulse with amplitude of roughly sech form, 
thus becoming reminiscent of NLS solitons. 
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II. THE DAVYDOV-SCOTT MODEL OF a-HELIX PROTEIN AND LATTICE 
NLS APPROXIMATION 

Davydov et al proposed a model of excitation propagation in a- helix protein in 2|. 
The version (171) of the model used here is based on refinements introduced by Scott in 3- 

n 

|5| which reflect newer experimental measurements. This model also neglects a variety of 
weaker interact terms described in [3j, and most physical parameters have been scaled away, 
leaving three time scales 

J ^ 1.4 THz, L ^ 2.3 THz, oj = ^Jxf^ ^ 12 THz. 



See also the surveys and references in , luj, |12|, and the similar earlier modeling by 
Holstein ^Q- 

The index n labels the amino acid residues. However, the helical structure is reflected by 
a repeating unit cell consisting of three consecutive residues, with total rotation of close to 
one full turn of the helix. Thus the physical structure is made clearer by indexing each unit 
cell by and the three residues within a unit cell with a = 0, 1, 2 (so that n = 3z/ + a), so 
occasionally the variables are labeled V'l^.o = and so on: 



dt 



The three substructures of residues with equal index a are almost linear and these are called 
spines: note that most interactions are along these spines. 



The Variables 

• The exciton variables ipn (coming from reduction of Schrodinger's equation) give the 
probability that the C=0 double bond in amino acid residue n is in a vibrational 
excited state. 

• The position variables g„ are the displacements of the residues from rest position, in 
the direction of the axis of the helix: that is, along spines. The corresponding momenta 
are Pn = mdqn/dt. 



• Parameter J measures the (attractive) interaction between excitons in residues that 
are adjacent along a spine. 

• Parameter L measures the (repulsive) interaction between excitons in residues that 
are adjacent along the helical molecular backbone. 

• Parameter m is the effective mass of each amino acid residue. Note that variations in 
mass between different amino acids are ignored: the oscillatory motion is assumed to 
be primarily in the part of the residue common to all amino acids, as the variable part 
is attached laterally at the outside of the helix. 

The Davydov-Scott System has a Hamiltonian form, with 

^ = Y1 -Ji'l'n^n+S + <+3^n) + ^(CV'n+l + <+lV^n) 
n 

+ Z2^+ 2^^"+' ~ + ~ ^nXV^n- (10) 

n 

Aside. The main difference in the original Davydov model was assuming a symmetrical 
form for the interaction between the exciton and mechanical variables: 

^(g„+3 - qn-3)lpn'^n- 
n 

Various boundary conditions can be imposed on the Hamiltonian. Defining the bond 
stretchings 

Sn '■= Qn+3 ~ Qn {= Qiy+l,a ~ Qu,a) , (H) 

zero end conditions are most natural physically: 

N residues indexed < n < N, with t/j^ — 0, s„ = for "out of bounds" values of the index 
n. 

For constructing PDE approximations via continuum limits, it is also convenient to con- 
sider an infinite chain with n e N and ^ 0, s„ ^ as \n\ oo. 

A. Approximation by a Lattice Nonlinear Schrodinger Equation 

Many phenomena of interest in solutions the Davydov system are retained in the approx- 
imation that the greater stiffness of the mechanical couplings causes the exciton quantities 
ipn to interact primarily with their time average, given by the singular limit m — > 0. 
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This gives the equation herein called the Davydov-Scott Lattice Nonlinear Schrddinger 
[DSLNLSj Equation 

+ + V-n-s) - L{i^n+1 + ^^n-l) + IV^nl'^n = 0. (12) 

This also has Hamiltonian form, with 

n n 

B. Momenta (Conserved Quantities Other Than The Hamiltonian) 

Charge. Both equations above have a conserved charge £ (also called exciton number 
or power depending on the physical application). This is related to the probability density 
of quantum mechanics, and notably, it is quadratic: 

n 

This charge is associated via Noether's Theorem with a linear symmetry group action, the 
gauge symmetry 

i^n^e^'i^n. rn^e-''rn- (15) 

Linear Momentum. The Davydov-Scott system also has a conserved momentum V, 
again (degenerately) quadratic: 

V = J2pn. (16) 

n 

This is associated via Noether's Theorem with the symmetry group action Qn ^ Qn + s. 
However, conservation of linear momentum is respected by almost any reasonable time 
discretization (for example, any Runge-Kutta method) so little more will be said about this. 

C. The NLS Continuum Limit 

The nonlinearity in ([8]) corresponds to an effective potential — I'j/'np, producing a "self- 
focusing" or "discrete self-trapping" effect that could counter dispersion and keep pulses 
compact and coherent over long propagation distances, akin to so-called solitons in PDE's. 
To seek this possible connection to solitary waves, a further approximation has been con- 
sidered by Davydov, Scott et al, based on the idea that solutions might evolve into a form 
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that is slowly varying along spines, due to the energetic favorability given by the intra-spine 
coupling terms, which can be rewritten as 



and \{qu+i,a - qu,aY- 

This will be seen below to require modification, but for comparison to previous work, 
note that the assumption of slowly varying form ipu,a{t) ~ ^/'^(z/Aa;, t) leads to a system of 
three linearly coupled nonlinear Schrodinger equations 

and then some further assumptions such as a phase shift form 

ipa = e'V, 30 = mod 27r 

and suitable rescaling leads to the focusing cubic NLS equation 

dip d^ip 
dt dx"^ 



with robust soliton solutions 



ip{x, t) = A sech(74(x — ft)) exp 



D. Momenta and Hamiltonian Symmetries via Invariant Quadratic Forms 

Conservation of momentum and charge can be verified directly (rather than invoking 
Noether's Theorem) by first noting that the state variables appear in the Hamiltonian only 
through the exciton products 

and the bond stretchings ( ITTi) . all of which are invariant under the associated gauge and 
translation symmetries, ensuring that the Hamiltonian has these symmetries. Then direct 
differentiation and application of the chain rule gives 

d£ dV 

— = 0, ^ = 0. 

dt ' dt 

However it will be convenient to use a more general calculation, given in section IIII B[ 



III. TIME DISCRETIZATION METHODS THAT CONSERVE ENERGY AND 
MOMENTA, AND RESPECT TIME-REVERSAL SYMMETRY 



The numerical methods will be described in terms of the more general Hamiltonian system 

'^^jvMy)-J^(y) (17) 

with J' an anti-symmetric matrix. One key to deriving methods that exactly conserve energy 
(the Hamiltonian) is to approximate (|T7|) with a discrete Hamiltonian system 

^ = JV,?/(y,y+) (18) 

by defining a suitable discrete gradient approximation 

Vy/(y,y+)^Vy/(y). (19) 

This can easily be done in ways that exactly conserve the Hamiltonian, but conserving 
other invariants (here all called momenta) requires an appropriate choice of the gradient 
approximation, and there is a natural limitation to quadratic (including linear) momenta. 



Some Notation for Difference Schemes 

We will focus on the time advance map for single time step, from a time t to t + h. 
Thus for a scalar variable y, a vector y, and likewise for other variables like q, p, and ip: 

• h = 6t denotes the change in t over the time step. 

• y alone without arguments denotes the value y{t) at time t, typically the beginning of 
the current time step. 

• t^ = t + h and y~^ denotes the value y{t^) = y{t + h). 

• Sy = y+-y. 

_ y + y^ 
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A. Energy Conserving discrete gradient methods 



The basic idea, originating in the work of 15|, |l6| is to define a discrete approximation 



Vy/ of the gradient depending on values at the beginning and end of a time step 

Vy/(y, y+) = (DJiy, y+), . . .) = J(y, y+), • • •) ^ Vy/(y), (20) 



and then solving the above discrete Hamiltonian system. In the case of a variable ma- 
trix i7(y), one instead uses an approximation i7(y,y^), with the natural choice being the 
midpoint approximation J7'(y,y^) = i7(y). We assume the natural consistency condition 

iim(Vy/)(y>y+) = Vy/(y). (21) 

A discrete gradient cannot simply be constructed from independently defined discrete 
approximations of the partial derivatives, because an important relation must be imposed 
on the components: all discrete gradients are required to satisfy the Discrete (Multivariable) 
Chain Rule 

'5/ = (Vy/)(y,y+)-5y. (22) 

This rule plus linearity will be imposed from now on. In fact, it will often be convenient to 
express discrete gradients through expressions for 6f in terms of (5y. 

Conservation of Energy With Discrete Gradient Methods 

Recall the proof of conservation of energy for a Hamiltonian system: dy/dt = J'Wy'H{y) 

— ^ = Vv'H • — multivariable chain rule 

dt ^ dt 

= Vy'H(y) • JVyl-Liy) Hamilton's equations 

= from the anti-symmetry of J . 

Conservation of energy is easily shown for a discrete gradient method by mimicking this 
argument: 

^ = (VyH)(y,y+) • ^ = (VyH)(y,y+) • J(VyH)(y,y+) = 0. 



B. Choosing a Discrete Gradient that Also Respects Quadratic Momenta 



On the other hand, conservation of momenta, such as the charge seen above, is not 
automatic. Instead one must choose amongst the infinity of possible discrete gradients to 
satisfy conservation laws. Here the approach introduced in [17|, ll8| is followed, based on 
three facts: 

1. There is a unique discrete gradient for functions of a single variable y 

Vyf{y,y^):=l (23) 
[^{y), y^ = y 

following from the chain rule requirement 

5f = ^yf{y,y^) 5y. 

2. There is a unique time-reversal symmetric discrete gradient for a product of two vari- 
ables 

SiVjVk) = VjSyk + Vk^Vj (24) 
which corresponds to evaluating the true gradient at the midpoint: 

V(%i/fc)(y,y+) = V(7/,?/fc)(y). (25) 

In fact this extends to a discrete product rule based on 

S{f9)=7S9 + 96f. (26) 

3. Many (perhaps all) physically relevant Hamiltonian systems with conserved momenta 
have a natural form in which all the momenta are quadratic (including linear) functions 
of the state variables, and are related through Noether's theorem to a group of affine 
symmetries of the Hamiltonian Ti, with invariance of H manifested by the fact that it 
can be expressed composition 

^(y) = ^(Q) (27) 

where each component Qi of the new state vector Q is a quadratic 

Qi = ^^^i'^yjyk + ^Kyj^ (each Ai ■= {Af } symmetric) (28) 



2 
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that is invariant under the symmetry group. For example, with the systems seen 
herein, the invariant quadratics with which the Hamiltonian can be expressed are the 
exciton products Hat and the bond stretchings s„. 

The discrete Jacobian of this change of variables is given by the true Jacobian evaluated 
at the midpoint: 

^yQ(y,y+) = ^yQ(y). 

These facts and the above chain rule requirement naturally lead to: 

A-^(y,y+) = 5^ AH(Q,Q+)AQ/(y'y+) = AH(Q,Q+)/^,gKy), (29) 

I I 

VyU = J^DiU VyQi, = Y,Di'H{Q, Q+) VyQi{y). (30) 

I I 

Using such a discrete gradient, energy and momentum will be conserved with any choice 
for the factors DiHiCl, Q"*")- In practice, the above rules for single variable functions, prod- 
ucts, compositions, and linearity are generally enough to construct a suitable discrete gra- 
dient for %. 



or 



Theorem 1 If for a Hamiltonian system f[T7j) . T-L has the form fl27j) as described above, 
and thus has a discrete gradient (!30|) . then numerical solution by the corresponding discrete 
gradient method ( JTSl) 

= JY. DiUiQ, Q+) V,Q, (^^) (31) 

conserves the Hamiltonian and all the quadratic momenta. 
The proof builds on Noether's Theorem via 



Lemma 1 Consider a Hamiltonian (1271) given in terms of quadratics fl28|) that are invariant 
under a continuous symmetry group, and let Q by any of the corresponding Noetherian 
invariants of the Hamiltonian flow, 

Q = \Y.A^\yk + Y.^yr (32) 

j,k j 

Then the following Poisson brackets vanish: 

{Q, Qi}{y) ■■= VQ(y) ■ JVQiiy) = 0. (33) 
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Proof of Lemma U\ The quadratic Q is invariant under the flow of any Hamiltonian 
H = H(Q), so in particular each choice Ti = Qi. gives a Hamiltonian flow with 

= ^ = VQ ■ JVn = VQ ■ JVQi = {Q, Qi}. (34) 

Proof of TheoremUl For any quadratic invariant of the original Hamiltonian flow, mim- 
icking ( l34l) gives 

^ = vg • jvn = VQ ■ jJ2(Din)VQi = 5^(ah(q, q+)){q, Qi}{y). (as) 

/ I 

Evaluating the Poisson brackets in (l33l) at y gives 6Q = 0. 



C. Practical Implementation: an Iterative Solution Method 

The system of equations will be nonlinear (unless the Hamiltonian system itself is linear), 
so we need an iterative solution method. The following method exploits the quasi-linearity 
of the system to preserves linear stability properties and exact momentum conservation 
without the cost of a full quasi-Newton method: Construct successive approximations y^'^^ 
of y"*" by solving 

^^^^ =jJ2 Di-HiQ, Q^'-'^) VyQi {^^^^ , (36) 

where Q*^'^-' := Q(y^'^''). Initialization can be with y*^*^) = y or some other suitable approx- 
imation of y^. That is, the nonlinear part Vq?^ is approximated using the current best 
available approximation y^'^"^) of y^, while the linear terms are left in terms of the un- 
known y'^'^^ to be solved for. This equation is linear in the unknown y'^^\ making its solution 
straightforward, and much as above, we have: 

Theorem 2 Each iterate y^'^^ given by the above scheme (l36ll conserves all quadratic first 
integrals that are conserved by the original discrete gradient scheme ( 13T]) . 

The proof is as for Theorem [T] above except that the Poisson brackets are evaluated at 

(y + y('')) /2. 

Energy is of course only conserved in the limit — oo. However, iterating until energy is 
accurate within machine rounding error is typically practical: if this take too many iterations, 
it is better for overall accuracy to reduce the time step size 6t to speed the convergence. 
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Unconditional Linear Stability Another advantage of this approach to iterative solution 
is that it has unconditional linear stability, since for a linear system, DqH is constant, the 
scheme converges in a single iteration, and is the A-stable implicit midpoint method. 



D. Higher Order Accuracy by Symmetric Step Composition 



The methods seen so far are only second order accurate in time. Fortunately, the method 
of symmetric step composition, developed in 19l-|22j for use with symplectic methods, and 
reviewed in the book 23[, gives a systematic way to construct methods of any higher even 
order while preserving all the interesting properties: conservation of the Hamiltonian and 
quadratic invariants, and time-reversal symmetry. 

Numerical results are computed below using the fourth-order accurate Suzuki method 



23 



, Example II. 4. 3, p. 45], composing five steps of lengths 7^5, 



7i = 72 = 74 = 75 



4- ^' 



73 = 1 - 471 



E. Aside: Comparison to Symplectic Methods 

The most commonly used conservative methods for Hamiltonian systems are symplectic 
methods such as the implicit midpoint method which conserve momenta, but cannot in 



general conserve energy, as described by a theorem of [2J]. For the systems here that are 
not in classical mechanical form 'H(q, p) = A'(p) + f/(q), the main choices are the midpoint 
method itself, high order diagonally implicit Runge-Kutta methods (which are all given by 
symmetric composition of midpoint steps), and fully implicit Gaussian methods. 

All but the last are cognates of the energy-momentum methods described here, reducing 
to the same form for linear equations, with the main difference being exact conservation of 
energy. Gaussian symplectic methods are desirable when the time step size is small enough 
to allow their solution by simple fixed point iteration, but are not cost effective for stiff 
systems where an unconditionally stable iterative method such as that above must be used. 
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IV. NUMERICAL RESULTS 



It seems most likely that the initial excitation will occur at a single residue so that the 
initial data is very far from the slowly varying form expected of a continuum limit. Thus the 
question addressed here, as in earlier work like , is whether such initial data is likely 
to lead to a form that can be well approximated by a smooth function of position, leading 
to a hopefully more tractable PDE model. It will be seen that some significant parts of the 
solution do fit a continuum limit, but instead of the integrable NLS model seen in the work 
of Davydov and Scott, they relate to a third derivative NLS equation and the Airy PDE, as 
discussed in Section |Vl 

A. Davydov-Scott System With Impulsive Initial Data 

The main case studied is initial data with excitation at a single residue uq: 



which seems physical plausible, given that the likely source of excitation is the transfer of 
energy from the ATP-ADP reaction to a vibrational excited state in one C=0 bond. Zero 
initial data is used for the mechanical terms: 



In addition to solving with the physically correct parameter values given above, the limit 
m — )■ is explored by comparing to solutions with far smaller value m = 10~^, and with the 
m = approximation DSLNLS in (|8]). 

B. Time step choice and handling of the faster mechanical terms 

The choice of time steps here is always constrained by 



which satisfies the natural accuracy and stability requirements for explicit methods and 
for convergence of a simple fixed point iterative method. While this is confirmed to give 
accurate solutions in the sense that all graphs here are completely indistinguishable from 



^„.(0) = A V-nlO) = forn ^ no. 



(37) 



gn(o) = p„(o) = 0. 
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results with any smaller time step size, far more is in fact true: with all graphs presented 
relating only the exciton variables, the graphical results are visually identical for any time 
step sizes 

6t < ^ 



2{J + L) 

depending only on the time scale manifested in the exciton evolution equation. Thus the 
time discretization is effectively handling any faster time scales in the mechanical variables 
in the innocuous way that one hopes for stiff modes to be handled by an unconditionally 
stable method, with no adverse effect on the accuracy of the more slowly evolving variables. 

C. Observations with an initial impulse at one residue 

Figures [1] and [2] show the results for initial data as in ( 137|) with A = 2 at the end of the 
molecule: n* = 0. (The "molecule" has an unnaturally large 3000 residues, to allow longer 
term trends to be seen more clearly.) 

A large portion of the initial excitation travels at roughly constant speed in a fairly 
coherent pulse, with far less dispersion than one would expect for a linear wave equation. 
In the forward half of the pulse, the amplitude is slowly varying, suggesting the possibility 
of a continuum limit description for part, but not all, of the solution. In these respects, the 
result is similar to those of js], which however used a 200 residue molecule, a shorter time, 
and the two-point initial impulse form ipoi^) = V'i(O) = A, ipniO) = for n > 1. Another 
noteworthy fact is that the region with slow change of amplitude between residues ends at 
almost exactly half way between the leading edge of the pulse and its origination point. This 
is as yet unexplained. 

However, slow variation is not seen in ipn, due to rapid phase variation: this is true even 
if one restricts to individual spines. This is seen by examining the real or imaginary parts 
along any spine or along all residues, as for example in Figure [3] for the spine a = 0. 

On closer examination, there is an approximate periodicity of period four in phase, which 
is revealed by examining the phase-shifted quantity Wn = i^ipn, as in Figure HI It is now 
truly slowly varying over the front half of the pulse, from residue 270 forward. This pattern 
is seen over a wide variety of impulsive initial data, varying amplitude A, using an impulse 
on the first one or two residues, and also with an initial impulse at another residue along 
the molecule. For example. Figure [5] shows a right-going pulse in w„ at t = 40, generated 
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by an initial impulse of amplitude A = 2 aX residue n* = 900. (There is a near identical 
left-going pulse, of course.) 



D. Smaller m and the limit m — t- 



The small mass and zero mass approximations work well. For example, with the same 
initial conditions as in Figures [THl] with the only change being setting m = 10~^ instead of 
its physical value 0.0069 (and reducing the time step size to 6t = 0.0005), the quantity Wn is 
seen in Figure El The values are almost identical to the corresponding ones for the m — )■ 
limit equation DSLNLS ([H]), as shown in Figure [3 In all cases: 

• Wherever there is slow variation in the amplitudes ipn, there is slow variation in residue 
index n of the Wn, and not of the ipn, even if one considers only values 4'm,a within a 
single spine a. 

• The behavior of w„ for the original Davydov-Scott system is well approximated by 
that in the LNLS approximation ([8]). 

• The leading pulse with slow variation moves at a speed of approximately 12.8 residues 
per unit time. 



V. THE THIRD DERIVATIVE NLS CONTINUUM LIMIT 

These observations can be partially explained by analysis of the dispersion relation for 
the linearization, and lead to consideration of a continuum limit for the transformation of 
([8]) into an equation for the w„. This is closely related to work of 25| and 26j on similar 



lattice NLS equations, with the main difference being that they address LNLS systems with 
only nearest-neighbor exciton dipole interactions: effectively J = 0. 

A. Dispersion Analysis and Phase Shift Along the Chain 



Consider the linear part of equation (IH 

.dip: 



dt 



+ J{lpn+3 + -ipn-s) - L{lpn+1 + 1pn-l) = 0. 



(38) 
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and seek discretized traveling waves of the form 



^„ = 0(2)e'(/^"+-*) , z = n-vt, 4>{z) 



with an emphasis on slow variation along the chain, indicated by small values of k. This 
gives the dispersion relation 

u;{k) = kv + 2Jcos{3{k + (3)) -2Lcos{k + (3). (39) 



As shown by Pelinovsky and Rothos 25|, using center manifold and normal form tech- 



niques (following earlier work of looss and Kirchgassner 27|, |28|] for the FPU and Klein- 
Gordon lattices) and using the "tail analysis" of Flach and Kladko 29(|, nonlinear solutions 
can bifurcate from solutions of this linear part at inflection points of ( 139|) . in particular 
uj'{0) = at /c = for the slowly varying solutions of interest here. This gives group velocity 

V = 6Jsin(3/3) -2Lsin(/3). 

Maximum signal speeds in each direction 

Vmax = ±(6J + 2L) [ = ±13 for the parameter values above] 

occur at the triple roots /3 = =F7r/2 of the dispersion relation. 

Aside: This is the main point where the detail of the lattice structure are important: 
specifically, that the "stride" in the coupling terms is odd. Behavior is different for a helix 
with an even number of residues per twist, though a near constant phase shift e^^ along the 
chain is still seen in numerical simulations, for somewhat different values of /3. 

The right-going case v = 6J + 2L, f3 = — 7r/2 corresponds to small amplitude solutions 
of the form tpn ~ (— i)"0(z)e^'^*. for which 

Wnit) = (i)"Vn(t) ~ e'(^(«~''*)+-*) (40) 

is slowly varying along the lattice. 



B. Approximation by the Third Derivative Nonlinear Schrodinger Equation 

Both the numerical observations and the above analysis suggest considering a continuum 
limit for the phase-shifted quantities Wn{t). Rewriting in terms of these gives 

-J-/ \t-/ XII*? / 1 \ 

+ J (w„+3 - Wn-z) + L[Wn+l - Wn-l) = l|Wn| (41) 
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as the starting point for such a continuum hmit. Next, expanding in the hmit of a slowly 
varying solution 

u!nit) ^ V^w{z,t), z = k{n — vt),T = et (42) 
where u = 6J + 2L as above, e = (1/3)kA;^, and k, = {27 J + L), gives 

Slow variation between consecutive nodes is indicated by A; ^ 1, which leads to approxima- 
tion by the third derivative nonlinear Schrddinger equation 

dw d^w , ,o , „. 

_ + _ . ,|.,|'^.,. (43) 

The continuum limit for the linearization considered above is 

dw d^w 

known as the Airy PDE (but also called the "Airy Diffusion Equation", though it not 
diffusive). This has a slowly decaying and spreading fundamental solution with initial 
impulse at z = 2:0, in terms of the |Airy function Ai] (see graph Airy_Functions.s"vg at 
IWikimedia Commonsjl : 



(3r)(i/3) V(3t)(i/3) 

Converting back via (H2|) with zq = kn^ to place the initial impulse at residue no gives 
approximate solution 

, , yk K^^^ . . f in — un) — vt\ , 
wJt) ^ V , Ai ^ — -^!f- . (46) 

C. Nonlinear Effects: Deviation from Airy Function Pulse Form and Evidence of 
Self-focusing 

The approximate solutions arising the above linearized continuum model have some no- 
table features: 

• Approximately traveling wave form, of speed f = 6 J + 2L = 13. 

• Airy function profile. 

• Slow decrease in amplitude and broadening, with time scale t^^^. 
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This speed is a good match for the values observed for the Davydov-Scott system, and 
the Airy function form is a good approximation for small initial amplitude, as is the slow 
temporal decrease in amplitude. 

The most noteworthy deviations as initial amplitude is increased are that the leading 
"hump" becomes more dominant (and so more similar to a scch pulse), and that its does 
not continue to spread in width and decrease in height at rate t^^^. Instead, it appears to 
stabilize in width. Perhaps this is the familiar control of dispersion by nonlinear self-focusing 
or "self-trapping" , but any such analysis remains to be done. 

VI. CONCLUSIONS 

1. The sustained traveling exciton pulses seen in Davydov-style exciton-oscillator model 
of energy propagation in a-hclix protein are well modeled by the "slow exciton" ap- 
proximation by a lattice NLS equation. 

2. As noted by other authors, the main part of the pulse has magnitude that varies 
slowly, allowing a PDE continuum limit approximation. 

3. However, the phase of the ipn varies rapidly in index n, by a quarter turn at each 
step, and thus the slow spatial variation is instead in Wn = (i)"'?/'n, leading to a new 
continuum limit approximation by the third derivative nonlinear Schrodinger equation 

4. Solutions for this bifurcate off from those of its linearization, the Airy PDE Wr+Wzzz = 
0, and retain a resemblance to the slowly modulated traveling wave solutions of that 
equation in terms of the Airy function even when the nonlinear term is of significant 
magnitude. 

5. Evidence of nonlinear "self- focusing" or "self-trapping" effects are seen, in that the 
leading hump of the pulse remains stronger and narrower as time increases than those 
of the linearization, supporting the more sustained propagation, despite the loss of the 
exact NLS sech soliton structure proposed in previous studies. 

6. The higher order exactly energy- momentum conserving time-discretization method 
used is seen to handle well the stiffness that can arise in such systems, making it a 
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good candidate for similar problems, including spatial discretization of various stiff 
nonlinear dispersive PDE's. 
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FIG. 1. Davydov-Scott system, A = 2 impulse initial data at residue 0: |V'n| at times t = 20,40. 
{St = 0.01.) 



0.25 




250 300 350 400 450 500 



residue index n 



FIG. 2. As above, t = 40 only: enlarged view of leading pulse amplitude. 
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FIG. 3. 3?(V'!/,o): real parts on spine a = 0; initial data as above. 
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FIG. 5. Real and imaginary parts of Wn = i"'V'n at t = 40 for initial impulse A = 2, n* = 900. 
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FIG. 6. As in Figure H] except with far smaller mass m = 10-^ 6t = 0.0005. 
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FIG. 7. As in Figures H] and [6] but for the DSLNLS equation. 
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